################################################################################
### Replication R script for Wellhausen & Zhu. "Exiting Russia." APSR.
####### R script to produce figures in the main text and appendix ##############

rm(list=ls(all=TRUE))

# install.packages("readstata13")
# install.packages("tidyverse")
# install.packages("ggpattern")
# install.packages("countrycode")
# install.packages("dotwhisker")
library(readstata13)
library(tidyverse)
library(ggpattern)
library(countrycode)
library(dotwhisker)

# set working directory
# setwd("")

df <- read.dta13("master_exit_russia.dta")
df <- df%>%filter(inactive==0&exit_removed==0&naics2022_2digit_text!="Public Administration")

################################################################################
######################### Figures in the Main Text #############################

## Distribution of Foreign-Invested Firm: Pre-Invasion Sample
# Figure 1(a): By Industry
df.ind2021 <- df%>%select(exit, naics2022_2digit_shorttext, naics2022_2digit_v2)%>%
  mutate(across(c(naics2022_2digit_shorttext), ~ na_if(., "")))%>%
  filter(is.na(naics2022_2digit_shorttext) == F&naics2022_2digit_v2!="92")%>%
  group_by(naics2022_2digit_shorttext, naics2022_2digit_v2)%>%
  count(naics2022_2digit_shorttext)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2021)%>%
  arrange(naics2022_2digit_v2)

df.ind2021 <- df.ind2021%>%mutate(naics = naics2022_2digit_shorttext,
                                  order = row_number())

ggplot(data=df.ind2021, aes(x= rev(factor(order)), y= count, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(unique(df.ind2021$naics)) ) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" )

ggsave(file="MNCs_byindustry2021.pdf",
           width=6,height=4)

# Figure 1(b): By Size: Pre-Invasion Sample
df.size2021 <- df%>%select(exit, size)%>%
  filter(is.na(size) == F)%>%
  group_by(size)%>%
  count(size)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2021)

ggplot(data=df.size2021, aes(x= rev(factor(size)), y= count, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 9.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(c("Small", "Medium", "Large", "Very Large") ) ) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" ) 

ggsave(file="MNCs_bysize2021.pdf",
       width=6,height=4)

rm(df.ind2021, df.size2021)

################################################################################
# Figure 2(a): Counts of Foreign-Invested Firms in Russia (2007--2021)
df2 <- read.dta13("exit_2007_2021.dta")%>%
  mutate(exit_pct = exit_rate*100,
         preinvasion = ifelse(year>=2021, 1, 0) )

ggplot(df2, aes(x=year, y = count_foreign,
               fill = as.factor(preinvasion),
               colour = as.factor(preinvasion)))+
  geom_bar(aes(y=count_foreign), stat="identity") +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(2,5,2,5), "cm")),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(margin = margin(t = 18, r = 5, b = -25, l = 5), size = 8, face = "bold", angle = 50),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(size = 8, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.5, 'cm'),
    legend.direction="vertical") +
  scale_x_continuous(breaks = c(2007:2021),labels= c(as.character(df2$year[1:14]), "Pre-Invasion" )) +
  labs(y = "Count",  x ="") 

ggsave(file="historical_counts.pdf",
    width=6,height=4.25)


# Figure 2(b): Exit Rates (2008--August 2023)
axis.x.labs <- c(paste0(df2$year[1:14],"–", c(as.character("08"), as.character("09"), as.character(c(10:21)))),
                 "Study Period")

ggplot(data=df2, aes(x= year, y=exit_pct)) +
  geom_point( aes(x= year, y=exit_pct, col = as.factor(preinvasion), shape = as.factor(preinvasion)),
              size = 3) +
  geom_line(aes(x= year, y=exit_pct, lty = as.factor(df2$preinvasion)),  color = "#F8766D") + 
  scale_linetype_manual(values=c(2, 0)) + 
  ylim(0,45) +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(2,5,2,5), "cm")),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(margin = margin(t = 18, r = 5, b = -25, l = 5), size = 8, face = "bold", angle = 50),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(size = 8, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.5, 'cm'),
    legend.direction="vertical")+
  scale_x_continuous(breaks = c(2007:2021),labels=  axis.x.labs) +
  labs(y = "Exit Rate in Percentage Points", x = "")

ggsave(file="historical_exits.pdf",
           width=6,height=4.25)

rm(df2)

################################################################################
##Figure 3: Subcompoents of the Exit Rate

df.exit <- df%>%select(exit_inactive, exit_f2d, exit_f2f)
exit.rate <- colMeans(df.exit, na.rm=TRUE)
exit.total <- colSums(df.exit, na.rm=TRUE)
exit.type <- c("Active to Inactive", "Foreign to Domestic", "Foreign to Foreign")

df.exit <- data.frame(exit.type, exit.rate, exit.total)

ggplot(data=df.exit, aes(x= exit.type, y=exit.rate*100, fill = exit.type)) +
  geom_bar(stat="identity") +
  geom_text(aes(label=exit.total, y = exit.rate*100 + .8), size = 4.5)+
  ylim(0,20) +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(2,5,0,5), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 0, l = 0), size = 8, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(size = 8, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="vertical")+
  labs(x = "", y = "Exit Rate in Percentage Points") 

ggsave(file="exit_type.pdf", width=4.5,height=4)

rm(df.exit, axis.x.labs, exit.rate, exit.total, exit.type)


################################################################################
# Figure 4: Exit Type

df.mlogit <- read.dta13("regcoef_mlogit.dta")%>%select(-N)%>%
  filter(substr(var, 1, 2)!="0:"&is.na(tstat)==F&grepl("naics2022", var)==F&
           substr(var,4, 7)!="cons")%>%mutate(equ = substr(var, 1, 1))%>%
  mutate(equ = recode(equ, "1" = "Active to Inactive",
                      "2" = "Foreign to Domestic",
                      "3" = "Foreign to Foreign") )%>%mutate(var = gsub("^(1|2|3):", "", var))

names(df.mlogit) <- c( "term", "estimate", "std.error", "statistic", "p.value", "model")

df.mlogit <- df.mlogit%>%filter(term%in%c("taxhaven", "eu", "uk","cyprus", "ukraine", "tradeflow_log", 
                                          "iia", "pta", "dtt", "2.size", "3.size", "4.size", "yrsinrus", "soe", 
                                          "listed_nonRU", "ru_subsidiaries", "mixed_nat", "multiple_homes", "multiple_fsh", "foreign_by_guo")==F)

dwplot(df.mlogit, vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2),
       dodge_size = .6)%>%
  relabel_predictors(c(unfriendly = "Unfriendly home state",
                       nslist = "On name-praise-shame list", 
                       consumer_oriented_3digit_median = "Consumer-oriented industry",
                       ruscontrol = "Russian managerial control",
                       russtrategic = "Russian strategic industry",
                       fs_stk_go_dum_median = "High fixed assets"
  ) )+
  theme_bw()+  
  theme( 
    plot.margin = margin(unit(c(2,5,5,5), "cm")),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text( size = 8, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=8, face = "bold"),
    legend.key.size = unit(.65, 'cm'),
    legend.direction="horizontal") + facet_grid(~model) +
  labs(x = "Estimates", y = "" ) 

ggsave(file="coefplot_mlogit.pdf", width=6,height=4)


################################################################################
######################### Figures in the Appendix ##############################

# Figure 1.1: Foreign-Invested Firms by "Unfriendly" Home States

df.unfriendly <- df%>%filter(unfriendly_bvd==1)%>%select(home)%>%
  group_by(home)%>% count(home, name = "frequency")%>%
  ungroup()%>%arrange(desc(frequency))%>%
  mutate(id = row_number())

ggplot(data=df.unfriendly, aes(x= rev(id), y= frequency, fill = "#F8766D")) +
  geom_bar(position="stack", stat="identity", width=.8) + 
  geom_text(aes(label=frequency), hjust= - .2, size = 2.8) + 
  coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(limits= rev(df.unfriendly$home)) +
  scale_y_continuous(limits = c(0, 9000), breaks = seq(from = 0, to = 9000, by = 3000)) +
  labs( y = "Frequency", x ="" )

ggsave(file="unfriendly_bystate.pdf", width=5,height=8)    

# Figure 1.2: Foreign-Invested Firms by "Friendly" Home States

df.friendly <- df%>%filter(unfriendly_bvd==0)%>%select(home)%>%
  group_by(home)%>% count(home, name = "frequency")%>%
  ungroup()%>%arrange(desc(frequency), home)%>%
  mutate(id = row_number())

ggplot(data=df.friendly, aes(x= rev(id), y= frequency, fill = "#F8766D")) +
  geom_bar(position="stack", stat="identity", width=.8) + 
  geom_text(aes(label=frequency), hjust= - .2, size = 2.8) + 
  coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 6, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(limits= rev(df.friendly$home)) +
  scale_y_continuous(limits = c(0, 2000), breaks = seq(from = 0, to = 2000, by = 500)) +
  labs( y = "Frequency", x ="" )

ggsave(file="friendly_bystate.pdf",  width=5,height=8)   

################################################################################
# Figure 3.4 “Name, Praise, and Shame” subsample comparison: By size

df.ns.size <- df%>%select(nslist, size)%>%filter(is.na(size) == F)%>%
  group_by(size)%>% count(nslist, name = "frequency")%>%
  ungroup()%>%arrange(nslist, size)%>%
  mutate(nslist2 = ifelse(nslist==0, "Not on list", "On name-praise-shame list"))

ggplot(data=df.ns.size, aes(x= size, y= frequency, fill = nslist2) ) +
  geom_bar(stat="identity", position = "dodge") + 
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(t= - 5, r = 5, b = -5, l = 5), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 10, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = c(.8, .92),
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.35, 'cm'),
    legend.direction="vertical") +
  scale_x_discrete(limits= c("Small", "Medium", "Large", "Very Large") ) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" )

ggsave(file="nslist_bysize.pdf", width=5,height=4)    

# Figure 3.5  “Name, Praise, and Shame” subsample comparison: By industry
df.ns.ind <- df%>%select(nslist, naics2022_2digit_shorttext, naics2022_2digit_v2)%>%
  mutate(across(c(naics2022_2digit_shorttext), ~ na_if(., "")))%>%
  filter(is.na(naics2022_2digit_shorttext) == F)%>%
  group_by(naics2022_2digit_v2, naics2022_2digit_shorttext)%>% count(nslist, name = "frequency")%>%
  arrange(naics2022_2digit_v2)%>% mutate(order = cur_group_id())%>%
  ungroup()%>%
  mutate(nslist2 = ifelse(nslist==0, "Not on list", "On name-praise-shame list"))

ggplot(data=df.ns.ind, aes(x= rev(factor(order)), y = (frequency), 
                           fill = (nslist2) ) )+
  geom_bar(stat="identity", position = "dodge") + coord_flip(ylim = c(0, 8000)) +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(t= - 5, r = 5, b = 5, l = 5), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 5, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = c(.75, .92),
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    #legend.spacing.y = unit(1, 'cm'),
    legend.direction="vertical") +
  scale_x_discrete(labels = rev(unique(df.ns.ind$naics2022_2digit_shorttext)) ) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" )

ggsave(file="nslist_byindustry.pdf", width=6,height=4)    


# Figure 3.6. “Name, Praise, and Shame” subsample comparison: By “unfriendly” home
df.ns.unfriendly <- df%>%filter(unfriendly_bvd==1)%>%select(nslist, home)%>%
  group_by(home, nslist)%>% count(nslist, name = "frequency")%>%
  ungroup()

df.ns.unfriendly <- left_join(df.ns.unfriendly, df.unfriendly, by = c("home"))%>%
  rename(frequency = frequency.x, total = frequency.y)%>%
  mutate(nslist2 = ifelse(nslist==0, "Not on list", "On name-praise-shame list"))%>%
  arrange(desc(id))

ggplot(data=df.ns.unfriendly, aes( x = fct_reorder(home, total), y= frequency, fill = nslist2) ) +
  geom_bar(position="dodge", stat="identity", width=.8) + 
  coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = c(.75, .75),
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="vertical") +
  labs(title="", x ="", y = "Count", colour = "", shape = "" )

ggsave(file="nslist_byhomeunfriendly.pdf", width=5,height=8)    


# Figure 3.7. “Name, Praise, and Shame” subsample comparison: By “friendly” home
df.ns.friendly <- df%>%filter(unfriendly_bvd==0)%>%select(nslist, home)%>%
  group_by(home, nslist)%>% count(nslist, name = "frequency")%>%
  ungroup()

df.ns.friendly <- left_join(df.ns.friendly, df.friendly, by = c("home"))%>%
  rename(frequency = frequency.x, total = frequency.y)%>%
  mutate(nslist2 = ifelse(nslist==0, "Not on list", "On mame-praise-shame list"))

ggplot(data=df.ns.friendly, aes( x = fct_reorder(home, total), y= frequency, fill = nslist2) ) +
  geom_bar(position="dodge", stat="identity", width=.8) + 
  coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 6, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = c(.70, .75),
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="vertical") +
  labs(title="", x ="", y = "Count", colour = "", shape = "" )

ggsave(file="nslist_byhomefriendly.pdf", width=5,height=8)    

rm(list = setdiff(ls(), "df"))

################################################################################
#############    Pre-invasion sample in historical context  #################
df.ind2021 <- df%>%select(exit, naics2022_2digit_shorttext, naics2022_2digit_v2)%>%
  mutate(across(c(naics2022_2digit_shorttext), ~ na_if(., "")))%>%
  filter(is.na(naics2022_2digit_shorttext) == F&naics2022_2digit_shorttext!="Public")%>%
  group_by(naics2022_2digit_shorttext, naics2022_2digit_v2)%>%
  count(naics2022_2digit_shorttext)%>%ungroup()%>%rename(count = n)%>%
  arrange(naics2022_2digit_v2)%>%mutate(year = 2021, order = row_number())

df.size2021 <- df%>%select(exit, size)%>%
  filter(is.na(size) == F)%>%
  group_by(size)%>%
  count(size)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2021)

df.home2021 <- df%>%select(exit, home_iso3)%>%filter(home_iso3!="")%>%
  mutate(region = countrycode(home_iso3, origin = 'iso3c', destination = 'region'))%>%
  group_by(region)%>% count(region)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2021)

df.naics2 <- df.ind2021%>%select(naics2022_2digit_shorttext, naics2022_2digit_v2)

# 2013--14

df.2013 <- read.dta13("placebo_exit_russia_2013-14.dta")

df.2013 <- df.2013%>%filter(inactive==0&naics2022_2digit_text!="Public Administration")

df.ind2013 <- df.2013%>%select(exit, naics2022_2digit_shorttext)%>%
  mutate(across(c(naics2022_2digit_shorttext), ~ na_if(., "")))%>%
  filter(is.na(naics2022_2digit_shorttext) == F&naics2022_2digit_shorttext!="Public")%>%
  group_by(naics2022_2digit_shorttext)%>%
  count(naics2022_2digit_shorttext)%>%ungroup()%>%rename(count = n)%>%
  mutate(year = 2013)

df.ind2013 <- left_join(df.ind2013, df.naics2, by = c("naics2022_2digit_shorttext"))%>%
  arrange(naics2022_2digit_v2)%>%mutate(order = row_number())

df.size2013 <- df.2013%>%select(exit, size)%>%
  filter(is.na(size) == F)%>%
  group_by(size)%>%
  count(size)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2013)

df.home2013 <- df.2013%>%select(exit, home_iso3)%>%filter(home_iso3!="")%>%
  mutate(region = countrycode(home_iso3, origin = 'iso3c', destination = 'region'))%>%
  group_by(region)%>% count(region)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2013)

# 2020--21

df.2020 <- read.dta13("placebo_exit_russia_2020-21.dta")

df.2020 <- df.2020%>%filter(inactive==0&naics2022_2digit_text!="Public Administration")

df.ind2020 <- df.2020%>%select(exit, naics2022_2digit_shorttext)%>%
  mutate(across(c(naics2022_2digit_shorttext), ~ na_if(., "")))%>%
  filter(is.na(naics2022_2digit_shorttext) == F&naics2022_2digit_shorttext!="Public")%>%
  group_by(naics2022_2digit_shorttext)%>%
  count(naics2022_2digit_shorttext)%>%ungroup()%>%rename(count = n)%>%
  mutate(year = 2020)

df.ind2020 <- left_join(df.ind2020, df.naics2, by = c("naics2022_2digit_shorttext"))%>%
  arrange(naics2022_2digit_v2)%>%mutate(order = row_number())

df.size2020 <- df.2020%>%select(exit, size)%>%
  filter(is.na(size) == F)%>%
  group_by(size)%>%
  count(size)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2020)

df.home2020 <- df.2020%>%select(exit, home_iso3)%>%filter(home_iso3!="")%>%
  mutate(region = countrycode(home_iso3, origin = 'iso3c', destination = 'region'))%>%
  group_by(region)%>% count(region)%>%ungroup()%>%rename(count = n)%>%mutate(year = 2020)

df.home <- rbind(df.home2021, df.home2020, df.home2013)

df.ind <- rbind(df.ind2021, df.ind2020, df.ind2013)%>%mutate(naics = naics2022_2digit_shorttext )

df.size <- rbind(df.size2021, df.size2020, df.size2013)

# Figure 4.8. Distribution of Foreign-Invested Firms by Industry (2013, 2020 & 2021)

ggplot(data=df.ind, aes(x= rev(factor(order)), y= count, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() + 
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 6.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 6, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(unique(df.ind$naics)) ) +
  scale_y_continuous(limits = c(0, 18000), breaks = seq(from = 0, to = 18000, by = 6000)) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" ) + facet_grid(~year)

ggsave(file="MNCs_byindustry.pdf", width=6,height=4)

# Figure 4.9. Distribution of Foreign-Invested Firms by Size (2013, 2020 & 2021)

ggplot(data=df.size, aes(x= rev(factor(size)), y= count, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 9.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(c("Small", "Medium", "Large", "Very Large") ) ) +
  scale_y_continuous(limits = c(0, 55009), breaks = seq(from = 0, to = 55009, by = 20000)) +
  labs(title="", x ="", y = "Count", colour = "", shape = "" ) + facet_grid(~year)

ggsave(file="MNCs_bysize.pdf", width=6,height=4)
  
# Figure 4.10. Distribution of Foreign-Invested Firms by Home Region (2013, 2020 & 2021)

ggplot(data=df.home, aes(x= rev(region), y= count, fill = "#F8766D")) +
  geom_bar( stat="identity", width = .8) + 
  coord_flip() +
  theme_bw()+ ylim(0, 53000)+
  theme(
    plot.margin = margin(unit(c(0,10,3,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 7.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 9.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.spacing.y = unit(1, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(df.home2020$region ) )+
  labs(title="", x ="", y = "Count", colour = "", shape = "" ) + facet_grid(~year)

ggsave(file="MNCs_byregion.pdf", width=6,height=4)

################################################################################
rm(list=ls(all=TRUE))

# Figure 8.11. Ownership data considerations, in historical comparison

df.miss <- read.dta13("owner_missingness.dta")

ggplot(data = df.miss, aes(x = year2)) +
  geom_point(aes(y = foreign_guo50, color = 'Foreign via GUO50 Pathway'), size = 3) +
  geom_point(aes(y = sh_miss, color = 'Shareholder Missingness'), size = 3, shape = 17) +
  geom_line(aes(y = foreign_guo50, color = 'Foreign via GUO50 Pathway'), linetype = 2) +
  geom_line(aes(y = sh_miss, color = 'Shareholder Missingness'), linetype = 2) +
  ylim(0, 40) +
  theme_bw() +
  theme(
    plot.margin = margin(unit(c(2, 5, 2, 5), "cm")),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(
      margin = margin(t = 18, r = 5, b = -23, l = 5),
      size = 8, face = "bold", angle = 50
    ),
    axis.title = element_text(size = 8, face = "bold"),
    axis.text.y = element_text(size = 8, face = "bold"),
    legend.position = c(0.2, 0.9),
    legend.text = element_text(size = 8, face = "bold"),
    legend.key.size = unit(0.6, 'cm'),
    legend.direction = "vertical"
  ) +
  scale_x_continuous(breaks = 1:16, labels = df.miss$year) +
  labs(y = "Percentage Points", x = "") +
  scale_color_manual(
    name = NULL,  
    breaks = c('Foreign via GUO50 Pathway', 'Shareholder Missingness'),
    values = c('Foreign via GUO50 Pathway' = '#F8766D', 'Shareholder Missingness' = '#00BFC4')
  ) + guides(
    color = guide_legend(
      override.aes = list(size = 2) 
    ) )

ggsave(file="sh_miss.pdf", width=6,height=4.25)

#-------------------------------------------------------------------------------
# Figure 8.12. Industry distribution of all firms in Russia (2021)
rm(list=ls(all=TRUE))

df <- read.dta13("firms_orbis_vs_fts_2021.dta")

ggplot(data= df, aes(x= rev(factor(order)), y= count_all/1000, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(0,2,2,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 6.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 6, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(unique(df$naics2022_2digit_shorttext)) ) +
  labs(title="", x ="", y = "Count (Thousands)", colour = "", shape = "" ) + facet_grid(~source)

ggsave(file="allfirms_byindustry2021.pdf", width=6,height=4)

# Figure 8.13. Industry distribution of foreign firms in Russia (2021)


ggplot(data= df, aes(x= rev(factor(order)), y= count_FIEs/1000, fill = "#F8766D")) +
  geom_bar(stat="identity") + coord_flip() +
  theme_bw()+
  theme(
    plot.margin = margin(unit(c(0,2,2,0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 6.5, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 0), size = 6, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="horizontal") +
  scale_x_discrete(labels= rev(unique(df$naics2022_2digit_shorttext)) ) +
  labs(title="", x ="", y = "Count (Thousands)", colour = "", shape = "" ) + facet_grid(~source)

ggsave(file="FIEs_byindustry2021.pdf", width=6,height=4)

rm(df)

################################################################################
# Figure 9.14. Percent of pre-invasion sample covered by economic treaties

df <- read.dta13("master_exit_russia.dta")%>%
  filter(inactive==0&exit_removed==0&naics2022_2digit_text!="Public Administration")

df.treaties <- df%>%select(pta, iia, dtt, unfriendly, home)%>%
  mutate(any = as.numeric(pta==1|iia==1|dtt==1),
         all = as.numeric(pta==1&iia==1&dtt==1))%>%group_by(unfriendly)%>%
  summarize(PTA = mean(pta, na.rm=TRUE),
            IIA =  mean(iia, na.rm=TRUE),
            DTT =  mean(dtt, na.rm=TRUE),
            Any =  mean(any, na.rm=TRUE),
            All =  mean(all, na.rm=TRUE))%>%ungroup()%>%
  pivot_longer(!unfriendly, names_to = "treaty", values_to = "pct")%>%
  mutate(home = ifelse(unfriendly==0, "Friendly", "Unfriendly"))

ggplot(data=df.treaties, aes( x = as.factor(rev(treaty)), y= pct*100, fill = treaty) ) +
  geom_bar(stat="identity", width=.9) +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(t= 3,r = 5,b=0,l=0), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = 2, l = 0), size = 8, face = "bold"),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0), size = 7.5, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=6, face = "bold"),
    legend.key.size = unit(.3, 'cm'),
    legend.direction="vertical") +
  scale_x_discrete(limits= c("PTA", "IIA", "DTT", "Any", "All") ) +
  facet_grid(~home) +
  labs(x ="", y = "" )

ggsave(file="Treaties_byunfriendly.pdf", width=6,height=4)  


################################################################################
# Figure 12.15: Determinants of foreign-invested firm exit: Comparing “unfriendly” and “friendly” home states

df.coef1 <- read.dta13("regcoef_unfriendly.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "Unfriendly")

df.coef2 <- read.dta13("regcoef_friendly.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "Friendly")

df.coef<-bind_rows(df.coef1, df.coef2)%>%
  select(var, coef, stderr, tstat, pval, model)%>%
  filter(var%in%c("mixed_nat", "multiple_homes", "multiple_fsh", "foreign_by_guo")==F)

names(df.coef) <- c("term", "estimate", "std.error", "statistic", "p.value", "model")

df.coef <- df.coef%>%filter(!term%in%c("eu", "ukraine", "uk", "pta", "2.size", "3.size", "4.size"))

dwplot(df.coef, vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2),
       dodge_size = .6, dot_args = list(aes(shape = model)),
       model_order = c("Unfriendly", "Friendly"))%>%
  relabel_predictors(c(nslist = "On name-praise-shame list", 
                       consumer_oriented_3digit_median = "Consumer-oriented industry",
                       ruscontrol = "Russian managerial control",
                       russtrategic = "Russian strategic industry",
                       fs_stk_go_dum_median = "High fixed assets",
                       tradeflow_log = "Trade with Russia",
                       taxhaven = "Tax haven",
                       cyprus = "Cyprus",
                       iia = "IIA", #pta = "PTA", 
                       dtt = "DTT",
                       yrsinrus = "Years in Russia",
                       ru_subsidiaries = "Russian subsidiaries",
                       soe = "State-owned",
                       listed_nonRU = "Listed outside of Russia") )+ 
  annotate(geom="text", x= -.68, y = 12, label="Hypotheses", angle = 90, size = 3)+
  annotate("segment",  x = -.64, y = 9.55, xend = -.64, yend = 14.5, size = .2) +
  annotate("segment",  x = -.64, y = 9.55, xend = -.63, yend = 9.55, size = .2) +
  annotate("segment",  x = -.64, y = 14.5, xend = -.63, yend = 14.5, size = .2) +
  annotate(geom="text", x= -.68, y = 7, label="Home Level",angle = 90, size = 3) +
  annotate("segment",  x = -.64, y = 4.55, xend = -.64, yend = 9.45, size = .2) +
  annotate("segment",  x = -.64, y = 9.45, xend = -.63, yend = 9.45, size = .2) +
  annotate("segment",  x = -.64, y = 4.55, xend = -.63, yend = 4.55, size = .2) +
  annotate(geom="text", x= -.68, y = 2.5, label="Firm Level",angle = 90, size = 3) +
  annotate("segment",  x = -.64, y = 0.5, xend = -.64, yend = 4.45, size = .2) +
  annotate("segment",  x = -.64, y = 4.45, xend = -.63, yend = 4.45, size = .2) +
  annotate("segment",  x = -.64, y = .5, xend = -.63, yend = .5, size = .2) +
  coord_cartesian(xlim = c(-.18, .26), clip = "off") +
  theme_bw()+ 
  theme(
    plot.margin = margin(unit(c(2,5,2,15), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = -20, l = 0), size = 8, face = "bold", angle = 0),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text(margin = margin(t = 2, r = 3, b = 5, l = 2), size = 8, face = "bold"),
    legend.position = "bottom",
    legend.title=element_blank(),
    legend.text=element_text(size=8, face = "bold"),
    legend.key.size = unit(.65, 'cm'),
    legend.direction="horizontal") + 
  scale_shape_manual(name = "Model", values = c(16, 17), breaks = c("Unfriendly", "Friendly")) +
  scale_colour_manual(name = "Model", values = c( "#F8766D", "#00BFC4") ,breaks = c("Unfriendly", "Friendly")) +
  guides( shape = guide_legend("Model"), colour = guide_legend("Model") ) +  
  geom_rect(data=data.frame(xmin=-Inf,xmax=Inf,ymin=-Inf,ymax= 4.5),
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
            fill = "deepskyblue1",alpha=0.1) +
  geom_rect(data=data.frame(xmin=-Inf,xmax=Inf,ymin=10.5,ymax= Inf),
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
            fill = "coral1",alpha=0.1)

ggsave(file="coefplot_unfriendly_vs_friendly.pdf",
    width=4,height=5)

rm(list=ls(all=TRUE))


################################################################################
# Figure 12.16: Determinants of foreign-invested  rm exit: Examining heterogeneity by size

df.coef.s <- read.dta13("regcoef_small.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "1")

df.coef.m <- read.dta13("regcoef_medium.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "2")

df.coef.l <- read.dta13("regcoef_large.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "3")

df.coef.xl <- read.dta13("regcoef_verylarge.dta")%>%
  filter(var!="1b.size"&var!="_cons")%>%mutate(model = "4")

df.coef.size<-bind_rows(df.coef.s, df.coef.m, df.coef.l, df.coef.xl)%>%
  select(var, coef, stderr, tstat, pval, model)%>%
  filter(var%in%c("mixed_nat", "multiple_homes", "multiple_fsh", "foreign_by_guo")==F)

names(df.coef.size) <- c("term", "estimate", "std.error", "statistic", "p.value", "model")

df.coef.size <- df.coef.size%>%filter(!term%in%c("eu", "ukraine", "uk", "pta", "2.size", "3.size", "4.size"))

###
annotation1 <- data.frame(
  model = "1", x= -1.45, y=12.5, label = "Hypotheses"
)

annotation2 <- data.frame(
  model = "1", x= -1.45, y = 7, label = "Home Level"
)

annotation3 <- data.frame(
  model = "1", x= -1.45, y = 2.5, label = "Firm Level"
)

segment1a <- data.frame(
  model = "1", x = -1.35, xend = -1.35, y = 9.55, yend = 15.5)
segment1b <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = 9.55, yend = 9.55)
segment1c <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = 15.5, yend = 15.5)

segment2a <- data.frame(
  model = "1", x = -1.35, xend = -1.35, y = 4.55, yend = 9.45)
segment2b <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = 4.55, yend = 4.55)
segment2c <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = 9.45, yend = 9.45)

segment3a <- data.frame(
  model = "1", x = -1.35, xend = -1.35, y = .5, yend = 4.45)
segment3b <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = .5, yend = .5)
segment3c <- data.frame(
  model = "1", x = -1.35, xend = -1.30, y = 4.45, yend = 4.45)

dwplot(df.coef.size, vline = geom_vline(xintercept = 0, colour = "grey50", linetype = 2),
       dodge_size = .6)%>%
  relabel_predictors(c(unfriendly = "Unfriendly home state",
                       nslist = "On name-praise-shame list", 
                       consumer_oriented_3digit_median = "Consumer-oriented industry",
                       ruscontrol = "Russian managerial control",
                       russtrategic = "Russian strategic industry",
                       fs_stk_go_dum_median = "High fixed assets",
                       tradeflow_log = "Trade with Russia",
                       taxhaven = "Tax haven",
                       cyprus = "Cyprus",
                       iia = "IIA", #pta = "PTA", 
                       dtt = "DTT",
                       yrsinrus = "Years in Russia",
                       ru_subsidiaries = "Russian subsidiaries",
                       soe = "State-owned",
                       listed_nonRU = "Listed outside of Russia") )+
  theme_bw()+
  facet_grid(~model, scales = "fixed", labeller = as_labeller(c('1' = "Small", '2' = "Medium", '3' = "Large", '4' = "Very Large"))) + 
  geom_text(data = annotation1, aes(x = x, y = y, label = label), size = 3, angle = 90) +
  geom_text(data = annotation2, aes(x = x, y = y, label = label), size = 3, angle = 90) +
  geom_text(data = annotation3, aes(x = x, y = y, label = label), size = 3, angle = 90) +
  geom_segment(data = segment1a,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment1b,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment1c,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment2a,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment2b,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment2c,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment3a,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment3b,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  geom_segment(data = segment3c,  aes(x = x, y = y, xend = xend, yend = yend), size = .2) +
  coord_cartesian(xlim = c(-.35, .3), clip = "off") +
  theme(
    plot.margin = margin(unit(c(2,5,2,15), "cm")),
    axis.text.x = element_text(margin = margin(t = 5, r = 0, b = -5, l = 0), size = 8, face = "bold", angle = 0),
    axis.title= element_text(size=8, face = "bold"),
    axis.text.y = element_text( size = 7, face = "bold"),
    legend.position = "none",
    legend.title=element_blank(),
    legend.text=element_text(size=8, face = "bold"),
    legend.key.size = unit(.65, 'cm'),
    legend.direction="horizontal") + 
    scale_x_continuous(breaks = seq(-.4, .2, .2), labels = seq(-.4, .2, .2)) +
   labs(x ="", y = "" ) +
  geom_rect(data=data.frame(xmin=-Inf,xmax=Inf,ymin=-Inf,ymax= 4.5),
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
            fill = "deepskyblue1",alpha=0.1) +
  geom_rect(data=data.frame(xmin=-Inf,xmax=Inf,ymin=9.5,ymax= Inf),
            aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
            fill = "coral1",alpha=0.1)

ggsave(file="coefplot_size.pdf",
    width=6.5,height=4)

